/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
**                                                                    **
** (C) Copyright 1999, The Regents of the University of California    **
** All Rights Reserved.                                               **
**                                                                    **
** Commercial use of this program without express permission of the   **
** University of California, Berkeley, is strictly prohibited.  See   **
** file 'COPYRIGHT'  in main directory for information on usage and   **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
**                                                                    **
** Developed by:                                                      **
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
**                                                                    **
** ****************************************************************** */

//
// 文件: GPUMatrixKernels.cu
// 创建时间: 2025年9月25日
// 说明: GPU矩阵组装CUDA核函数实现
// 作者: OpenSees GPU加速项目组
//

// ============================================================
// Performance optimization: Disable kernel debug output
// ============================================================
#ifndef GPU_ENABLE_DEBUG_OUTPUT
#define GPU_ENABLE_DEBUG_OUTPUT 0
#endif

#if GPU_ENABLE_DEBUG_OUTPUT
#define KERNEL_DEBUG_PRINT(...) KERNEL_DEBUG_PRINT(__VA_ARGS__)
#else
#define KERNEL_DEBUG_PRINT(...) ((void)0)
#endif

#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cuda_runtime_api.h>
#include <cublas_v2.h>
#include <cusolverDn.h>
#include <math.h>
#include <stdio.h>

#include "GPUElementSupport.h"

// CUDA错误检查宏
#define CUDA_CHECK(call) \
    do { \
        cudaError_t error = call; \
        if (error != cudaSuccess) { \
            KERNEL_DEBUG_PRINT("CUDA Error in %s at line %d: %s\n", __FILE__, __LINE__, cudaGetErrorString(error)); \
            return -1; \
        } \
    } while(0)

// 数学常量定义
__device__ const double PI = 3.14159265358979323846;
__device__ const double SQRT_2 = 1.41421356237309504880;
__device__ const double SQRT_3 = 1.73205080756887729353;

// Double precision atomicAdd for older compute capabilities
// 注意：CUDA 6.0+ (compute capability 6.0+) 原生支持double的atomicAdd
// 对于更早的架构，需要使用atomicCAS实现
#if defined(__CUDA_ARCH__) && __CUDA_ARCH__ < 600
// For older compute capabilities (< 6.0), implement atomicAdd for double using atomicCAS
__device__ double atomicAdd(double* address, double val) {
    unsigned long long int* address_as_ull = (unsigned long long int*)address;
    unsigned long long int old = *address_as_ull, assumed;

    do {
        assumed = old;
        old = atomicCAS(address_as_ull, assumed,
                        __double_as_longlong(val + __longlong_as_double(assumed)));
    } while (assumed != old);

    return __longlong_as_double(old);
}
#endif
// For compute capability 6.0 and higher, native atomicAdd for double is available from CUDA

// GPU全局内存中的共享数据指针声明（在GPUSolver.cpp中分配）
extern __device__ double* d_elementStiffness;
extern __device__ double* d_elementResidual;
extern __device__ int* d_elementDOF;
extern __device__ int* d_elementTypes;
extern __device__ double* d_elementData;

//==============================================================================
// CUDA核函数：计算Truss单元刚度矩阵
//==============================================================================

__global__ void computeTrussStiffness_kernel(
    int numElements,                  // 单元数量
    const GPUElementData* elementData,// 单元数据数组
    double* stiffnessMatrix,         // 输出：单元刚度矩阵 [numElements × maxDOF²]
    double* tempMatrix               // 临时计算空间
) {
    int elementIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (elementIdx >= numElements) return;
    
    // 获取当前单元数据
    const GPUElementData& element = elementData[elementIdx];
    
    // 检查是否为Truss单元
    GPUSupportedElementType elemType = static_cast<GPUSupportedElementType>(element.elementType);
    if (elemType != GPUSupportedElementType::TRUSS && 
        elemType != GPUSupportedElementType::TRUSS_SECTION &&
        elemType != GPUSupportedElementType::TRUSS2) {
        return;
    }
    
    // 获取材料参数
    double E = element.materialProps[0];    // 弹性模量
    double A = element.materialProps[1];    // 截面积
    // double rho = element.materialProps[2];  // 密度（质量矩阵时使用）
    
    // 获取节点坐标
    double x1 = element.nodeCoords[0], y1 = element.nodeCoords[1], z1 = element.nodeCoords[2];
    double x2 = element.nodeCoords[3], y2 = element.nodeCoords[4], z2 = element.nodeCoords[5];
    
    // 计算单元长度
    double dx = x2 - x1;
    double dy = y2 - y1;
    double dz = z2 - z1;
    double L = sqrt(dx*dx + dy*dy + dz*dz);
    
    if (L <= 1e-12) return;  // 避免长度为零的单元
    
    // 计算方向余弦
    double cosX = dx / L;
    double cosY = dy / L;
    double cosZ = dz / L;
    
    // 计算轴向刚度
    double k_axial = E * A / L;
    
    // 构建局部坐标系刚度矩阵（轴向刚度）
    // k_local = k_axial * [1  -1]
    //                     [-1  1]
    
    // 构建坐标变换矩阵 T (6x6)
    // T = [cosX cosY cosZ   0     0     0  ]
    //     [  0    0    0  cosX  cosY  cosZ]
    
    // 全局坐标系刚度矩阵 K = T^T * k_local * T
    // 计算结果：K是6x6矩阵，元素为：
    // K_ij = k_axial * (T[0,i] * T[0,j] - T[0,i] * T[1,j] - T[1,i] * T[0,j] + T[1,i] * T[1,j])
    
    // 获取输出矩阵位置（假设maxDOFPerElement = 72）
    const int maxDOFPerElement = 72;
    double* K = &stiffnessMatrix[elementIdx * maxDOFPerElement * maxDOFPerElement];
    
    // 初始化刚度矩阵为0
    for (int i = 0; i < GPUElementConstants::TRUSS_STIFFNESS_SIZE; i++) {
        K[i] = 0.0;
    }
    
    // 填充6×6刚度矩阵（3D Truss单元）
    // DOF排列：[u1, v1, w1, u2, v2, w2]
    
    // 节点1的贡献 (DOF 0,1,2)
    K[0*6 + 0] = k_axial * cosX * cosX;     // K11 = k * cosX^2
    K[0*6 + 1] = k_axial * cosX * cosY;     // K12 = k * cosX * cosY
    K[0*6 + 2] = k_axial * cosX * cosZ;     // K13 = k * cosX * cosZ
    K[0*6 + 3] = -k_axial * cosX * cosX;    // K14 = -k * cosX^2
    K[0*6 + 4] = -k_axial * cosX * cosY;    // K15 = -k * cosX * cosY
    K[0*6 + 5] = -k_axial * cosX * cosZ;    // K16 = -k * cosX * cosZ
    
    K[1*6 + 0] = k_axial * cosY * cosX;     // K21 = k * cosY * cosX
    K[1*6 + 1] = k_axial * cosY * cosY;     // K22 = k * cosY^2
    K[1*6 + 2] = k_axial * cosY * cosZ;     // K23 = k * cosY * cosZ
    K[1*6 + 3] = -k_axial * cosY * cosX;    // K24 = -k * cosY * cosX
    K[1*6 + 4] = -k_axial * cosY * cosY;    // K25 = -k * cosY^2
    K[1*6 + 5] = -k_axial * cosY * cosZ;    // K26 = -k * cosY * cosZ
    
    K[2*6 + 0] = k_axial * cosZ * cosX;     // K31 = k * cosZ * cosX
    K[2*6 + 1] = k_axial * cosZ * cosY;     // K32 = k * cosZ * cosY
    K[2*6 + 2] = k_axial * cosZ * cosZ;     // K33 = k * cosZ^2
    K[2*6 + 3] = -k_axial * cosZ * cosX;    // K34 = -k * cosZ * cosX
    K[2*6 + 4] = -k_axial * cosZ * cosY;    // K35 = -k * cosZ * cosY
    K[2*6 + 5] = -k_axial * cosZ * cosZ;    // K36 = -k * cosZ^2
    
    // 节点2的贡献 (DOF 3,4,5)
    K[3*6 + 0] = -k_axial * cosX * cosX;    // K41 = -k * cosX^2
    K[3*6 + 1] = -k_axial * cosX * cosY;    // K42 = -k * cosX * cosY
    K[3*6 + 2] = -k_axial * cosX * cosZ;    // K43 = -k * cosX * cosZ
    K[3*6 + 3] = k_axial * cosX * cosX;     // K44 = k * cosX^2
    K[3*6 + 4] = k_axial * cosX * cosY;     // K45 = k * cosX * cosY
    K[3*6 + 5] = k_axial * cosX * cosZ;     // K46 = k * cosX * cosZ
    
    K[4*6 + 0] = -k_axial * cosY * cosX;    // K51 = -k * cosY * cosX
    K[4*6 + 1] = -k_axial * cosY * cosY;    // K52 = -k * cosY^2
    K[4*6 + 2] = -k_axial * cosY * cosZ;    // K53 = -k * cosY * cosZ
    K[4*6 + 3] = k_axial * cosY * cosX;     // K54 = k * cosY * cosX
    K[4*6 + 4] = k_axial * cosY * cosY;     // K55 = k * cosY^2
    K[4*6 + 5] = k_axial * cosY * cosZ;     // K56 = k * cosY * cosZ
    
    K[5*6 + 0] = -k_axial * cosZ * cosX;    // K61 = -k * cosZ * cosX
    K[5*6 + 1] = -k_axial * cosZ * cosY;    // K62 = -k * cosZ * cosY
    K[5*6 + 2] = -k_axial * cosZ * cosZ;    // K63 = -k * cosZ^2
    K[5*6 + 3] = k_axial * cosZ * cosX;     // K64 = k * cosZ * cosX
    K[5*6 + 4] = k_axial * cosZ * cosY;     // K65 = k * cosZ * cosY
    K[5*6 + 5] = k_axial * cosZ * cosZ;     // K66 = k * cosZ^2
}

//==============================================================================
// CUDA核函数：计算梁单元刚度矩阵
//==============================================================================

__global__ void computeBeamStiffness_kernel(
    int numElements,                  // 单元数量
    const GPUElementData* elementData,// 单元数据数组
    double* stiffnessMatrix,         // 输出：单元刚度矩阵
    double* tempMatrix               // 临时计算空间
) {
    int elementIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (elementIdx >= numElements) return;
    
    // 获取当前单元数据
    const GPUElementData& element = elementData[elementIdx];
    
    // 检查是否为梁单元
    GPUSupportedElementType elemType = static_cast<GPUSupportedElementType>(element.elementType);
    bool is2DBeam = (elemType == GPUSupportedElementType::ELASTIC_BEAM_2D ||
                     elemType == GPUSupportedElementType::MOD_ELASTIC_BEAM_2D ||
                     elemType == GPUSupportedElementType::ELASTIC_TIMOSHENKO_BEAM_2D);
                     
    bool is3DBeam = (elemType == GPUSupportedElementType::ELASTIC_BEAM_3D ||
                     elemType == GPUSupportedElementType::MOD_ELASTIC_BEAM_3D ||
                     elemType == GPUSupportedElementType::ELASTIC_TIMOSHENKO_BEAM_3D);
    
    if (!is2DBeam && !is3DBeam) return;
    
    // 获取材料参数
    double E  = element.materialProps[0];    // 弹性模量
    double G  = element.materialProps[1];    // 剪切模量
    double A  = element.materialProps[2];    // 截面积
    double Iy = element.materialProps[3];    // y轴惯性矩
    double Iz = element.materialProps[4];    // z轴惯性矩
    double J  = element.materialProps[5];    // 扭转惯性矩
    // double rho = element.materialProps[6]; // 密度
    // double alpha = element.materialProps[7]; // 热膨胀系数
    
    // 获取节点坐标
    double x1 = element.nodeCoords[0], y1 = element.nodeCoords[1], z1 = element.nodeCoords[2];
    double x2 = element.nodeCoords[3], y2 = element.nodeCoords[4], z2 = element.nodeCoords[5];
    
    // 计算单元长度
    double dx = x2 - x1;
    double dy = y2 - y1;
    double dz = z2 - z1;
    double L = sqrt(dx*dx + dy*dy + dz*dz);
    
    if (L <= 1e-12) return;  // 避免长度为零的单元
    
    // 获取输出矩阵位置
    const int maxDOFPerElement = 72;
    double* K = &stiffnessMatrix[elementIdx * maxDOFPerElement * maxDOFPerElement];
    
    if (is2DBeam) {
        // 2D梁单元刚度矩阵 (6×6)
        // DOF排列：[u1, v1, theta1, u2, v2, theta2]
        
        // 初始化为0
        for (int i = 0; i < GPUElementConstants::BEAM_2D_STIFFNESS_SIZE; i++) {
            K[i] = 0.0;
        }
        
        // 计算刚度参数
        double EA_L = E * A / L;                    // 轴向刚度
        double EI_L = E * Iy / L;                   // 弯曲刚度基数
        double EI_L2 = EI_L / L;                    // 弯曲刚度/L
        double EI_L3 = EI_L2 / L;                   // 弯曲刚度/L²
        
        // 方向余弦
        double cosTheta = dx / L;
        double sinTheta = dy / L;
        
        // 局部坐标系刚度矩阵
        double k_local[36] = {0.0}; // 6×6矩阵
        
        // 轴向刚度项
        k_local[0*6 + 0] = EA_L;        // K11 = EA/L
        k_local[0*6 + 3] = -EA_L;       // K14 = -EA/L
        k_local[3*6 + 0] = -EA_L;       // K41 = -EA/L
        k_local[3*6 + 3] = EA_L;        // K44 = EA/L
        
        // 弯曲刚度项
        k_local[1*6 + 1] = 12.0 * EI_L3;    // K22 = 12EI/L³
        k_local[1*6 + 2] = 6.0 * EI_L2;     // K23 = 6EI/L²
        k_local[1*6 + 4] = -12.0 * EI_L3;   // K25 = -12EI/L³
        k_local[1*6 + 5] = 6.0 * EI_L2;     // K26 = 6EI/L²
        
        k_local[2*6 + 1] = 6.0 * EI_L2;     // K32 = 6EI/L²
        k_local[2*6 + 2] = 4.0 * EI_L;      // K33 = 4EI/L
        k_local[2*6 + 4] = -6.0 * EI_L2;    // K35 = -6EI/L²
        k_local[2*6 + 5] = 2.0 * EI_L;      // K36 = 2EI/L
        
        k_local[4*6 + 1] = -12.0 * EI_L3;   // K52 = -12EI/L³
        k_local[4*6 + 2] = -6.0 * EI_L2;    // K53 = -6EI/L²
        k_local[4*6 + 4] = 12.0 * EI_L3;    // K55 = 12EI/L³
        k_local[4*6 + 5] = -6.0 * EI_L2;    // K56 = -6EI/L²
        
        k_local[5*6 + 1] = 6.0 * EI_L2;     // K62 = 6EI/L²
        k_local[5*6 + 2] = 2.0 * EI_L;      // K63 = 2EI/L
        k_local[5*6 + 4] = -6.0 * EI_L2;    // K65 = -6EI/L²
        k_local[5*6 + 5] = 4.0 * EI_L;      // K66 = 4EI/L
        
        // 坐标变换矩阵 T (6×6)
        double T[36] = {0.0};
        T[0*6 + 0] = cosTheta;  T[0*6 + 1] = sinTheta;   // T11, T12
        T[1*6 + 0] = -sinTheta; T[1*6 + 1] = cosTheta;   // T21, T22
        T[2*6 + 2] = 1.0;                                // T33
        T[3*6 + 3] = cosTheta;  T[3*6 + 4] = sinTheta;   // T44, T45
        T[4*6 + 3] = -sinTheta; T[4*6 + 4] = cosTheta;   // T54, T55
        T[5*6 + 5] = 1.0;                                // T66
        
        // 计算全局刚度矩阵 K = T^T * k_local * T
        // 这里使用简化计算，直接给出结果
        for (int i = 0; i < 6; i++) {
            for (int j = 0; j < 6; j++) {
                double sum = 0.0;
                for (int k = 0; k < 6; k++) {
                    for (int l = 0; l < 6; l++) {
                        sum += T[k*6 + i] * k_local[k*6 + l] * T[l*6 + j];
                    }
                }
                K[i*6 + j] = sum;
            }
        }
    }
    else if (is3DBeam) {
        // 3D梁单元刚度矩阵 (12×12) - 这是一个更复杂的实现
        // DOF排列：[u1, v1, w1, rx1, ry1, rz1, u2, v2, w2, rx2, ry2, rz2]
        
        // 初始化为0
        for (int i = 0; i < GPUElementConstants::BEAM_3D_STIFFNESS_SIZE; i++) {
            K[i] = 0.0;
        }
        
        // 简化实现：只考虑轴向刚度，弯曲和扭转刚度需要更复杂的坐标变换
        double EA_L = E * A / L;
        double cosX = dx / L, cosY = dy / L, cosZ = dz / L;
        
        // 轴向刚度项 (简化版本)
        K[0*12 + 0] = EA_L * cosX * cosX;
        K[0*12 + 6] = -EA_L * cosX * cosX;
        K[6*12 + 0] = -EA_L * cosX * cosX;
        K[6*12 + 6] = EA_L * cosX * cosX;
        
        // TODO: 完整的3D梁刚度矩阵实现（包括弯曲和扭转）
        // 这需要复杂的坐标变换和局部坐标系处理
    }
}

//==============================================================================
// CUDA核函数：计算四边形单元刚度矩阵
//==============================================================================

__global__ void computeQuadStiffness_kernel(
    int numElements,                  // 单元数量
    const GPUElementData* elementData,// 单元数据数组
    double* stiffnessMatrix,         // 输出：单元刚度矩阵
    double* tempMatrix               // 临时计算空间
) {
    int elementIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (elementIdx >= numElements) return;
    
    // 获取当前单元数据
    const GPUElementData& element = elementData[elementIdx];
    
    // 检查是否为四边形单元
    GPUSupportedElementType elemType = static_cast<GPUSupportedElementType>(element.elementType);
    if (elemType != GPUSupportedElementType::FOUR_NODE_QUAD && 
        elemType != GPUSupportedElementType::FOUR_NODE_QUAD_3D &&
        elemType != GPUSupportedElementType::FOUR_NODE_QUAD_02) {
        return;
    }
    
    // 获取材料参数
    double E = element.materialProps[0];      // 弹性模量
    double nu = element.materialProps[1];     // 泊松比
    double thickness = element.materialProps[2]; // 厚度
    // double rho = element.materialProps[3];  // 密度
    
    // 获取节点坐标（4个节点）
    double x1 = element.nodeCoords[0], y1 = element.nodeCoords[1];
    double x2 = element.nodeCoords[3], y2 = element.nodeCoords[4]; 
    double x3 = element.nodeCoords[6], y3 = element.nodeCoords[7];
    double x4 = element.nodeCoords[9], y4 = element.nodeCoords[10];
    
    // 获取输出矩阵位置
    const int maxDOFPerElement = 72;
    double* K = &stiffnessMatrix[elementIdx * maxDOFPerElement * maxDOFPerElement];
    
    // 初始化刚度矩阵为0
    for (int i = 0; i < GPUElementConstants::QUAD_STIFFNESS_SIZE; i++) {
        K[i] = 0.0;
    }
    
    // 平面应力弹性矩阵
    double factor = E / (1.0 - nu * nu);
    double D[9]; // 3×3弹性矩阵
    D[0] = factor;                  // D11 = E/(1-ν²)
    D[1] = factor * nu;             // D12 = νE/(1-ν²)
    D[2] = 0.0;                     // D13 = 0
    D[3] = factor * nu;             // D21 = νE/(1-ν²)
    D[4] = factor;                  // D22 = E/(1-ν²)
    D[5] = 0.0;                     // D23 = 0
    D[6] = 0.0;                     // D31 = 0
    D[7] = 0.0;                     // D32 = 0
    D[8] = factor * (1.0 - nu) / 2.0; // D33 = G = E/[2(1+ν)]
    
    // 高斯积分点（2×2积分）
    const double gaussPts[2] = {-0.577350269189626, 0.577350269189626};  // ±1/√3
    const double gaussWts[2] = {1.0, 1.0};
    
    // 对每个积分点进行积分
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            double xi = gaussPts[i];
            double eta = gaussPts[j];
            double wt = gaussWts[i] * gaussWts[j];
            
            // 形函数及其导数
            double N[4], dN_dxi[4], dN_deta[4];
            
            // 双线性形函数
            N[0] = 0.25 * (1.0 - xi) * (1.0 - eta);    // N1
            N[1] = 0.25 * (1.0 + xi) * (1.0 - eta);    // N2
            N[2] = 0.25 * (1.0 + xi) * (1.0 + eta);    // N3
            N[3] = 0.25 * (1.0 - xi) * (1.0 + eta);    // N4
            
            // 形函数对自然坐标的导数
            dN_dxi[0] = -0.25 * (1.0 - eta);           // dN1/dξ
            dN_dxi[1] =  0.25 * (1.0 - eta);           // dN2/dξ
            dN_dxi[2] =  0.25 * (1.0 + eta);           // dN3/dξ
            dN_dxi[3] = -0.25 * (1.0 + eta);           // dN4/dξ
            
            dN_deta[0] = -0.25 * (1.0 - xi);           // dN1/dη
            dN_deta[1] = -0.25 * (1.0 + xi);           // dN2/dη
            dN_deta[2] =  0.25 * (1.0 + xi);           // dN3/dη
            dN_deta[3] =  0.25 * (1.0 - xi);           // dN4/dη
            
            // 计算Jacobian矩阵
            double J11 = dN_dxi[0]*x1 + dN_dxi[1]*x2 + dN_dxi[2]*x3 + dN_dxi[3]*x4;
            double J12 = dN_dxi[0]*y1 + dN_dxi[1]*y2 + dN_dxi[2]*y3 + dN_dxi[3]*y4;
            double J21 = dN_deta[0]*x1 + dN_deta[1]*x2 + dN_deta[2]*x3 + dN_deta[3]*x4;
            double J22 = dN_deta[0]*y1 + dN_deta[1]*y2 + dN_deta[2]*y3 + dN_deta[3]*y4;
            
            // Jacobian行列式
            double detJ = J11 * J22 - J12 * J21;
            if (fabs(detJ) < 1e-12) continue; // 避免奇异单元
            
            // Jacobian逆矩阵
            double invJ11 =  J22 / detJ;
            double invJ12 = -J12 / detJ;
            double invJ21 = -J21 / detJ;
            double invJ22 =  J11 / detJ;
            
            // 形函数对全局坐标的导数
            double dN_dx[4], dN_dy[4];
            for (int k = 0; k < 4; k++) {
                dN_dx[k] = invJ11 * dN_dxi[k] + invJ12 * dN_deta[k];
                dN_dy[k] = invJ21 * dN_dxi[k] + invJ22 * dN_deta[k];
            }
            
            // 构建应变-位移矩阵B (3×8)
            // B = [dN1/dx    0    dN2/dx    0    dN3/dx    0    dN4/dx    0  ]
            //     [   0   dN1/dy    0   dN2/dy    0   dN3/dy    0   dN4/dy]
            //     [dN1/dy dN1/dx dN2/dy dN2/dx dN3/dy dN3/dx dN4/dy dN4/dx]
            
            double B[24]; // 3×8矩阵
            for (int k = 0; k < 4; k++) {
                B[0*8 + 2*k]     = dN_dx[k];  // εxx项
                B[0*8 + 2*k + 1] = 0.0;
                B[1*8 + 2*k]     = 0.0;       // εyy项
                B[1*8 + 2*k + 1] = dN_dy[k];
                B[2*8 + 2*k]     = dN_dy[k];  // γxy项
                B[2*8 + 2*k + 1] = dN_dx[k];
            }
            
            // 计算B^T * D * B并累加到刚度矩阵
            double volume = thickness * detJ * wt;  // 体积元素
            
            for (int m = 0; m < 8; m++) {
                for (int n = 0; n < 8; n++) {
                    double sum = 0.0;
                    for (int p = 0; p < 3; p++) {
                        for (int q = 0; q < 3; q++) {
                            sum += B[p*8 + m] * D[p*3 + q] * B[q*8 + n];
                        }
                    }
                    K[m*8 + n] += volume * sum;
                }
            }
        }
    }
}

//==============================================================================
// CUDA核函数：全局刚度矩阵组装
//==============================================================================

__global__ void assembleStiffnessMatrix_kernel(
    int numElements,                  // 单元数量
    const GPUElementData* elementData,// 单元数据数组
    const double* elementStiffness,   // 单元刚度矩阵
    double* globalStiffness,          // 全局刚度矩阵（压缩存储格式）
    int globalSize                    // 全局系统大小
) {
    int elementIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (elementIdx >= numElements) return;
    
    // 获取当前单元数据
    const GPUElementData& element = elementData[elementIdx];
    
    // 获取单元DOF数量
    int elementDOF = element.numDOF;
    if (elementDOF <= 0) return;
    
    // 获取单元刚度矩阵
    const int maxDOFPerElement = 72;
    const double* K_elem = &elementStiffness[elementIdx * maxDOFPerElement * maxDOFPerElement];
    
    // 组装到全局矩阵（原子操作确保线程安全）
    for (int i = 0; i < elementDOF; i++) {
        int globalI = element.DOFIndices[i];
        if (globalI < 0 || globalI >= globalSize) continue;
        
        for (int j = 0; j < elementDOF; j++) {
            int globalJ = element.DOFIndices[j];
            if (globalJ < 0 || globalJ >= globalSize) continue;
            
            // 计算全局矩阵位置（假设为满矩阵存储）
            double value = K_elem[i * maxDOFPerElement + j];
            if (fabs(value) > 1e-16) {
                atomicAdd(&globalStiffness[globalI * globalSize + globalJ], value);
            }
        }
    }
}

//==============================================================================
// CUDA核函数：计算Truss单元残差向量
//==============================================================================

__global__ void computeTrussResidual_kernel(
    int numElements,                  // 单元数量
    const GPUElementData* elementData,// 单元数据数组
    const double* displacement,       // 节点位移向量
    double* residualVector,          // 输出：单元残差向量 [numElements × maxDOF]
    double* tempVector               // 临时计算空间
) {
    int elementIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (elementIdx >= numElements) return;
    
    // 获取当前单元数据
    const GPUElementData& element = elementData[elementIdx];
    
    // 检查是否为Truss单元
    GPUSupportedElementType elemType = static_cast<GPUSupportedElementType>(element.elementType);
    if (elemType != GPUSupportedElementType::TRUSS && 
        elemType != GPUSupportedElementType::TRUSS_SECTION &&
        elemType != GPUSupportedElementType::TRUSS2) {
        return;
    }
    
    // 获取材料参数
    double E = element.materialProps[0];    // 弹性模量
    double A = element.materialProps[1];    // 截面积
    // 当前应力由CPU传入，存储在stateVariables中
    double stress = element.stateVariables[0];   // 当前轴向应力（CPU计算后传入）
    
    // 获取节点坐标
    double x1 = element.nodeCoords[0], y1 = element.nodeCoords[1], z1 = element.nodeCoords[2];
    double x2 = element.nodeCoords[3], y2 = element.nodeCoords[4], z2 = element.nodeCoords[5];
    
    // 计算单元长度和方向余弦
    double dx = x2 - x1;
    double dy = y2 - y1;
    double dz = z2 - z1;
    double L = sqrt(dx*dx + dy*dy + dz*dz);
    
    if (L <= 1e-12) return;  // 避免长度为零的单元
    
    double cosX = dx / L;
    double cosY = dy / L;
    double cosZ = dz / L;
    
    // 获取单元DOF位移
    const int maxDOFPerElement = 72;
    double u_elem[6]; // Truss单元6个DOF [u1, v1, w1, u2, v2, w2]
    
    for (int i = 0; i < 6; i++) {
        int globalDOF = element.DOFIndices[i];
        if (globalDOF >= 0) {
            u_elem[i] = displacement[globalDOF];
        } else {
            u_elem[i] = 0.0;
        }
    }
    
    // 计算轴向位移（局部坐标系）
    double u_axial = (u_elem[3] - u_elem[0]) * cosX + 
                     (u_elem[4] - u_elem[1]) * cosY + 
                     (u_elem[5] - u_elem[2]) * cosZ;
    
    // 计算轴向应变
    double strain = u_axial / L;
    
    // 计算轴向力（如果没有材料非线性，使用线性关系）
    double axialForce;
    if (fabs(stress) > 1e-16) {
        // 使用CPU传入的应力
        axialForce = stress * A;
    } else {
        // 使用线性弹性关系
        axialForce = E * A * strain;
    }
    
    // 计算节点等效力（全局坐标系）
    // 节点1的力（与轴向力方向相反）
    double F1x = -axialForce * cosX;
    double F1y = -axialForce * cosY;
    double F1z = -axialForce * cosZ;
    
    // 节点2的力（与轴向力方向相同）
    double F2x = axialForce * cosX;
    double F2y = axialForce * cosY;
    double F2z = axialForce * cosZ;
    
    // 获取输出残差向量位置
    double* R = &residualVector[elementIdx * maxDOFPerElement];
    
    // 初始化残差向量为0
    for (int i = 0; i < 6; i++) {
        R[i] = 0.0;
    }
    
    // 填充残差向量
    R[0] = F1x;  // 节点1 X方向力
    R[1] = F1y;  // 节点1 Y方向力
    R[2] = F1z;  // 节点1 Z方向力
    R[3] = F2x;  // 节点2 X方向力
    R[4] = F2y;  // 节点2 Y方向力
    R[5] = F2z;  // 节点2 Z方向力
}

//==============================================================================
// CUDA核函数：计算梁单元残差向量
//==============================================================================

__global__ void computeBeamResidual_kernel(
    int numElements,                  // 单元数量
    const GPUElementData* elementData,// 单元数据数组
    const double* displacement,       // 节点位移向量
    double* residualVector,          // 输出：单元残差向量
    double* tempVector               // 临时计算空间
) {
    int elementIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (elementIdx >= numElements) return;
    
    // 获取当前单元数据
    const GPUElementData& element = elementData[elementIdx];
    
    // 检查是否为梁单元
    GPUSupportedElementType elemType = static_cast<GPUSupportedElementType>(element.elementType);
    bool is2DBeam = (elemType == GPUSupportedElementType::ELASTIC_BEAM_2D ||
                     elemType == GPUSupportedElementType::MOD_ELASTIC_BEAM_2D ||
                     elemType == GPUSupportedElementType::ELASTIC_TIMOSHENKO_BEAM_2D);
                     
    bool is3DBeam = (elemType == GPUSupportedElementType::ELASTIC_BEAM_3D ||
                     elemType == GPUSupportedElementType::MOD_ELASTIC_BEAM_3D ||
                     elemType == GPUSupportedElementType::ELASTIC_TIMOSHENKO_BEAM_3D);
    
    if (!is2DBeam && !is3DBeam) return;
    
    // 获取材料参数
    double E  = element.materialProps[0];    // 弹性模量
    double A  = element.materialProps[2];    // 截面积
    double Iy = element.materialProps[3];    // y轴惯性矩
    
    // 获取当前内力（CPU计算后传入）
    double axialForce = element.stateVariables[0];     // 轴向力
    double shearForce = element.stateVariables[1];     // 剪切力
    double moment1 = element.stateVariables[2];        // 节点1弯矩
    double moment2 = element.stateVariables[3];        // 节点2弯矩
    
    // 获取节点坐标
    double x1 = element.nodeCoords[0], y1 = element.nodeCoords[1];
    double x2 = element.nodeCoords[3], y2 = element.nodeCoords[4];
    
    // 计算单元长度和方向
    double dx = x2 - x1;
    double dy = y2 - y1;
    double L = sqrt(dx*dx + dy*dy);
    
    if (L <= 1e-12) return;
    
    double cosTheta = dx / L;
    double sinTheta = dy / L;
    
    // 获取输出残差向量位置
    const int maxDOFPerElement = 72;
    double* R = &residualVector[elementIdx * maxDOFPerElement];
    
    if (is2DBeam) {
        // 2D梁单元残差向量 (6个DOF)
        // DOF排列：[u1, v1, theta1, u2, v2, theta2]
        
        // 初始化残差向量
        for (int i = 0; i < 6; i++) {
            R[i] = 0.0;
        }
        
        // 如果没有CPU传入的内力，使用简化线性计算
        if (fabs(axialForce) <= 1e-16 && fabs(shearForce) <= 1e-16) {
            // 获取单元位移
            double u_elem[6];
            for (int i = 0; i < 6; i++) {
                int globalDOF = element.DOFIndices[i];
                if (globalDOF >= 0) {
                    u_elem[i] = displacement[globalDOF];
                } else {
                    u_elem[i] = 0.0;
                }
            }
            
            // 局部坐标系位移
            double u1_local = cosTheta * u_elem[0] + sinTheta * u_elem[1];
            double v1_local = -sinTheta * u_elem[0] + cosTheta * u_elem[1];
            double theta1_local = u_elem[2];
            double u2_local = cosTheta * u_elem[3] + sinTheta * u_elem[4];
            double v2_local = -sinTheta * u_elem[3] + cosTheta * u_elem[4];
            double theta2_local = u_elem[5];
            
            // 轴向应变和力
            double axial_strain = (u2_local - u1_local) / L;
            axialForce = E * A * axial_strain;
            
            // 弯曲变形和内力（简化计算）
            double EI = E * Iy;
            double d_v = v2_local - v1_local;
            double d_theta = theta2_local - theta1_local;
            
            shearForce = 12.0 * EI / (L * L * L) * d_v + 6.0 * EI / (L * L) * (theta1_local + theta2_local);
            moment1 = 4.0 * EI / L * theta1_local + 2.0 * EI / L * theta2_local - 6.0 * EI / (L * L) * d_v;
            moment2 = 2.0 * EI / L * theta1_local + 4.0 * EI / L * theta2_local - 6.0 * EI / (L * L) * d_v;
        }
        
        // 局部坐标系节点力
        double F1_local[3] = {-axialForce, -shearForce, moment1};
        double F2_local[3] = {axialForce, shearForce, moment2};
        
        // 转换到全局坐标系
        R[0] = cosTheta * F1_local[0] - sinTheta * F1_local[1];  // F1x
        R[1] = sinTheta * F1_local[0] + cosTheta * F1_local[1];  // F1y
        R[2] = F1_local[2];                                       // M1
        R[3] = cosTheta * F2_local[0] - sinTheta * F2_local[1];  // F2x
        R[4] = sinTheta * F2_local[0] + cosTheta * F2_local[1];  // F2y
        R[5] = F2_local[2];                                       // M2
    }
    else if (is3DBeam) {
        // 3D梁单元残差向量 (12个DOF) - 简化实现
        // DOF排列：[u1, v1, w1, rx1, ry1, rz1, u2, v2, w2, rx2, ry2, rz2]
        
        // 初始化残差向量
        for (int i = 0; i < 12; i++) {
            R[i] = 0.0;
        }
        
        // 简化为轴向力计算（弯曲和扭转待完善）
        if (fabs(axialForce) <= 1e-16) {
            double u_elem[12];
            for (int i = 0; i < 12; i++) {
                int globalDOF = element.DOFIndices[i];
                if (globalDOF >= 0) {
                    u_elem[i] = displacement[globalDOF];
                } else {
                    u_elem[i] = 0.0;
                }
            }
            
            double dz = element.nodeCoords[5] - element.nodeCoords[2];
            double L3D = sqrt(dx*dx + dy*dy + dz*dz);
            double cosZ = dz / L3D;
            
            // 轴向位移
            double u_axial = (u_elem[6] - u_elem[0]) * cosTheta + 
                             (u_elem[7] - u_elem[1]) * sinTheta + 
                             (u_elem[8] - u_elem[2]) * cosZ;
            
            axialForce = E * A * u_axial / L3D;
        }
        
        // 轴向力分量
        R[0] = -axialForce * cosTheta;   // F1x
        R[1] = -axialForce * sinTheta;   // F1y
        R[6] = axialForce * cosTheta;    // F2x
        R[7] = axialForce * sinTheta;    // F2y
        
        // TODO: 完整的3D梁内力计算（弯曲和扭转）
    }
}

//==============================================================================
// CUDA核函数：计算四边形单元残差向量
//==============================================================================

__global__ void computeQuadResidual_kernel(
    int numElements,                  // 单元数量
    const GPUElementData* elementData,// 单元数据数组
    const double* displacement,       // 节点位移向量
    double* residualVector,          // 输出：单元残差向量
    double* tempVector               // 临时计算空间
) {
    int elementIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (elementIdx >= numElements) return;
    
    // 获取当前单元数据
    const GPUElementData& element = elementData[elementIdx];
    
    // 检查是否为四边形单元
    GPUSupportedElementType elemType = static_cast<GPUSupportedElementType>(element.elementType);
    if (elemType != GPUSupportedElementType::FOUR_NODE_QUAD && 
        elemType != GPUSupportedElementType::FOUR_NODE_QUAD_3D &&
        elemType != GPUSupportedElementType::FOUR_NODE_QUAD_02) {
        return;
    }
    
    // 获取材料参数
    double E = element.materialProps[0];      // 弹性模量
    double nu = element.materialProps[1];     // 泊松比
    double thickness = element.materialProps[2]; // 厚度
    
    // 获取当前应力状态（CPU计算后传入）
    double sigma_xx = element.stateVariables[0];   // 正应力 σxx
    double sigma_yy = element.stateVariables[1];   // 正应力 σyy
    double tau_xy = element.stateVariables[2];     // 剪应力 τxy
    
    // 获取节点坐标
    double x1 = element.nodeCoords[0], y1 = element.nodeCoords[1];
    double x2 = element.nodeCoords[3], y2 = element.nodeCoords[4]; 
    double x3 = element.nodeCoords[6], y3 = element.nodeCoords[7];
    double x4 = element.nodeCoords[9], y4 = element.nodeCoords[10];
    
    // 获取单元位移
    double u_elem[8]; // 4个节点，每个节点2个DOF
    for (int i = 0; i < 8; i++) {
        int globalDOF = element.DOFIndices[i];
        if (globalDOF >= 0) {
            u_elem[i] = displacement[globalDOF];
        } else {
            u_elem[i] = 0.0;
        }
    }
    
    // 获取输出残差向量位置
    const int maxDOFPerElement = 72;
    double* R = &residualVector[elementIdx * maxDOFPerElement];
    
    // 初始化残差向量
    for (int i = 0; i < 8; i++) {
        R[i] = 0.0;
    }
    
    // 高斯积分点（2×2积分）
    const double gaussPts[2] = {-0.577350269189626, 0.577350269189626};
    const double gaussWts[2] = {1.0, 1.0};
    
    // 如果没有CPU传入的应力，计算线性弹性应力
    bool useLinearStress = (fabs(sigma_xx) + fabs(sigma_yy) + fabs(tau_xy) <= 1e-16);
    
    // 平面应力弹性矩阵
    double factor = E / (1.0 - nu * nu);
    double D[9];
    D[0] = factor;                          // D11
    D[1] = factor * nu;                     // D12
    D[2] = 0.0;                             // D13
    D[3] = factor * nu;                     // D21
    D[4] = factor;                          // D22
    D[5] = 0.0;                             // D23
    D[6] = 0.0;                             // D31
    D[7] = 0.0;                             // D32
    D[8] = factor * (1.0 - nu) / 2.0;      // D33
    
    // 对每个积分点计算内力
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            double xi = gaussPts[i];
            double eta = gaussPts[j];
            double wt = gaussWts[i] * gaussWts[j];
            
            // 形函数导数
            double dN_dxi[4], dN_deta[4];
            dN_dxi[0] = -0.25 * (1.0 - eta);
            dN_dxi[1] =  0.25 * (1.0 - eta);
            dN_dxi[2] =  0.25 * (1.0 + eta);
            dN_dxi[3] = -0.25 * (1.0 + eta);
            
            dN_deta[0] = -0.25 * (1.0 - xi);
            dN_deta[1] = -0.25 * (1.0 + xi);
            dN_deta[2] =  0.25 * (1.0 + xi);
            dN_deta[3] =  0.25 * (1.0 - xi);
            
            // Jacobian矩阵
            double J11 = dN_dxi[0]*x1 + dN_dxi[1]*x2 + dN_dxi[2]*x3 + dN_dxi[3]*x4;
            double J12 = dN_dxi[0]*y1 + dN_dxi[1]*y2 + dN_dxi[2]*y3 + dN_dxi[3]*y4;
            double J21 = dN_deta[0]*x1 + dN_deta[1]*x2 + dN_deta[2]*x3 + dN_deta[3]*x4;
            double J22 = dN_deta[0]*y1 + dN_deta[1]*y2 + dN_deta[2]*y3 + dN_deta[3]*y4;
            
            double detJ = J11 * J22 - J12 * J21;
            if (fabs(detJ) < 1e-12) continue;
            
            // Jacobian逆矩阵
            double invJ11 =  J22 / detJ;
            double invJ12 = -J12 / detJ;
            double invJ21 = -J21 / detJ;
            double invJ22 =  J11 / detJ;
            
            // 形函数对全局坐标的导数
            double dN_dx[4], dN_dy[4];
            for (int k = 0; k < 4; k++) {
                dN_dx[k] = invJ11 * dN_dxi[k] + invJ12 * dN_deta[k];
                dN_dy[k] = invJ21 * dN_dxi[k] + invJ22 * dN_deta[k];
            }
            
            // 应变-位移矩阵B (3×8)
            double B[24];
            for (int k = 0; k < 4; k++) {
                B[0*8 + 2*k]     = dN_dx[k];
                B[0*8 + 2*k + 1] = 0.0;
                B[1*8 + 2*k]     = 0.0;
                B[1*8 + 2*k + 1] = dN_dy[k];
                B[2*8 + 2*k]     = dN_dy[k];
                B[2*8 + 2*k + 1] = dN_dx[k];
            }
            
            // 计算当前积分点的应力
            double stress[3];
            if (useLinearStress) {
                // 计算应变
                double strain[3] = {0.0, 0.0, 0.0};
                for (int k = 0; k < 8; k++) {
                    strain[0] += B[0*8 + k] * u_elem[k];  // εxx
                    strain[1] += B[1*8 + k] * u_elem[k];  // εyy
                    strain[2] += B[2*8 + k] * u_elem[k];  // γxy
                }
                
                // 计算应力 σ = D * ε
                stress[0] = D[0] * strain[0] + D[1] * strain[1];  // σxx
                stress[1] = D[3] * strain[0] + D[4] * strain[1];  // σyy
                stress[2] = D[8] * strain[2];                      // τxy
            } else {
                // 使用CPU传入的应力
                stress[0] = sigma_xx;
                stress[1] = sigma_yy;
                stress[2] = tau_xy;
            }
            
            // 计算节点等效力 F = B^T * σ * volume
            double volume = thickness * detJ * wt;
            for (int k = 0; k < 8; k++) {
                double force = 0.0;
                for (int l = 0; l < 3; l++) {
                    force += B[l*8 + k] * stress[l];
                }
                R[k] += volume * force;
            }
        }
    }
}

//==============================================================================
// CUDA核函数：全局残差向量组装
//==============================================================================

__global__ void assembleResidualVector_kernel(
    int numElements,                  // 单元数量
    const GPUElementData* elementData,// 单元数据数组
    const double* elementResidual,    // 单元残差向量
    double* globalResidual,           // 全局残差向量
    int globalSize                    // 全局系统大小
) {
    int elementIdx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (elementIdx >= numElements) return;
    
    // 获取当前单元数据
    const GPUElementData& element = elementData[elementIdx];
    
    // 获取单元DOF数量
    int elementDOF = element.numDOF;
    if (elementDOF <= 0) return;
    
    // 获取单元残差向量
    const int maxDOFPerElement = 72;
    const double* R_elem = &elementResidual[elementIdx * maxDOFPerElement];
    
    // 组装到全局向量（原子操作确保线程安全）
    for (int i = 0; i < elementDOF; i++) {
        int globalI = element.DOFIndices[i];
        if (globalI < 0 || globalI >= globalSize) continue;
        
        double value = R_elem[i];
        if (fabs(value) > 1e-16) {
            atomicAdd(&globalResidual[globalI], value);
        }
    }
}

//==============================================================================
// C++封装函数：供GPUSolver调用
//==============================================================================

extern "C" {

//==============================================================================
// 残差向量计算核函数的C++封装调用函数
// 注意：刚度矩阵计算的启动函数在GPUElementKernels.cu中定义
//==============================================================================

// 调用Truss单元残差计算核函数
int launchTrussResidualKernel(
    int numElements,
    const GPUElementData* d_elementData,
    const double* d_displacement,
    double* d_residualVector,
    double* d_tempVector,
    cudaStream_t stream = 0
) {
    if (numElements <= 0) return 0;
    
    // 计算网格配置
    int blockSize = 256;
    int gridSize = (numElements + blockSize - 1) / blockSize;
    
    // 启动核函数
    computeTrussResidual_kernel<<<gridSize, blockSize, 0, stream>>>(
        numElements, d_elementData, d_displacement, d_residualVector, d_tempVector
    );
    
    // 检查CUDA错误
    cudaError_t error = cudaGetLastError();
    if (error != cudaSuccess) {
        KERNEL_DEBUG_PRINT("CUDA Truss Residual Kernel Error: %s\n", cudaGetErrorString(error));
        return -1;
    }
    
    return 0;
}

// 调用梁单元残差计算核函数
int launchBeamResidualKernel(
    int numElements,
    const GPUElementData* d_elementData,
    const double* d_displacement,
    double* d_residualVector,
    double* d_tempVector,
    cudaStream_t stream = 0
) {
    if (numElements <= 0) return 0;
    
    int blockSize = 256;
    int gridSize = (numElements + blockSize - 1) / blockSize;
    
    computeBeamResidual_kernel<<<gridSize, blockSize, 0, stream>>>(
        numElements, d_elementData, d_displacement, d_residualVector, d_tempVector
    );
    
    cudaError_t error = cudaGetLastError();
    if (error != cudaSuccess) {
        KERNEL_DEBUG_PRINT("CUDA Beam Residual Kernel Error: %s\n", cudaGetErrorString(error));
        return -1;
    }
    
    return 0;
}

// 调用四边形单元残差计算核函数
int launchQuadResidualKernel(
    int numElements,
    const GPUElementData* d_elementData,
    const double* d_displacement,
    double* d_residualVector,
    double* d_tempVector,
    cudaStream_t stream = 0
) {
    if (numElements <= 0) return 0;
    
    int blockSize = 256;
    int gridSize = (numElements + blockSize - 1) / blockSize;
    
    computeQuadResidual_kernel<<<gridSize, blockSize, 0, stream>>>(
        numElements, d_elementData, d_displacement, d_residualVector, d_tempVector
    );
    
    cudaError_t error = cudaGetLastError();
    if (error != cudaSuccess) {
        KERNEL_DEBUG_PRINT("CUDA Quad Residual Kernel Error: %s\n", cudaGetErrorString(error));
        return -1;
    }
    
    return 0;
}

// 调用全局残差向量组装核函数
int launchResidualAssemblyKernel(
    int numElements,
    const GPUElementData* d_elementData,
    const double* d_elementResidual,
    double* d_globalResidual,
    int globalSize,
    cudaStream_t stream = 0
) {
    if (numElements <= 0 || globalSize <= 0) return 0;
    
    int blockSize = 256;
    int gridSize = (numElements + blockSize - 1) / blockSize;
    
    assembleResidualVector_kernel<<<gridSize, blockSize, 0, stream>>>(
        numElements, d_elementData, d_elementResidual, d_globalResidual, globalSize
    );
    
    cudaError_t error = cudaGetLastError();
    if (error != cudaSuccess) {
        KERNEL_DEBUG_PRINT("CUDA Residual Assembly Kernel Error: %s\n", cudaGetErrorString(error));
        return -1;
    }
    
    return 0;
}

} // extern "C"